function [M,PI,G,ylist,xlist,lpd]=resolkw(system,track,comp,param,gainkeep);


answinv = param(9);
whatinfo = param(10);

% Rational expectations solution program implementing 
% the various steps described in R.G.King and M.W.Watson
% "System Reduction and Model Solution Algorithms for
% Linear Rational Expectations Models." The inputs and 
% outputs are summarized by the function statement:
%
%     [M,PI,G,ylist,xlist,lpd]=resolkw(system,driver,track,comp)
%
% There are two required inputs to the program:
%
%   (1) The file 'system.m' must describe the 
%       linear rational expectations model to be solved.
%       (the input file need not be called 'system.m', but
%        can take any name that the user supplies in the call
%        to the resolkw.m program).  The RE model takes the form:
%
%         A Ey(t+1)|I(t) = B y(t) + C0 x(t) + ... Cn Ex(t+n)|I(t)
%
%       system.m must contain 
%         (a) the matrices A,B,CF=[C0 ... Cn]
%         (b) a vector of predetermined variable locations (lpd)
%         (c) the number of exogenous variables, 
%              which implies the number of leads (nlead=cols(CF)/nx)
%         (d) lists of endogenous (ylist) and exogenous (xlist) variables
%
%   (2) The file 'driver.m' must contain the matrices RHO and Q
%       that describe the driving processes for the system:
%         x(t) = Q delta(t)
%         delta(t) = RHO * delta(t-1) + Gi e(t)
%       The matrix G of this driving process is not essential for
%       the current solution.
% 
%   (3) The program output is a solution in state space form:
%       
%           y(t) = PI * s(t)
%          
%           s(t) = M s(t-1) + G e(t)
%      
%       with s(t) containing the predetermined endogenous variables
%       followed by the delta(t).
%
%   (4) The user may optionally decide not to track the solution of the
%       model (except for forced error messages). This involves setting 
%       track=0 in the call to redkw.m
%
%   (5) The user may optionally decide not to use the Schur decomposition
%       method.  This involves supply a two element vector that involves
%       [schur nobalance].  The parameter schur = 0 involves computation via
%       the eigenvector-eigenvalue method; the parameter nobalance = 1
%       involves the computation of the eigenvalues by MATLAB's no balance
%       option, which leads to more accurate computations in some situations.

if (nargin<2)
   track=1;   % tracking is default
end

% --------------------------------------------------------------------------
% Select details of computational method: 
%    standard eigenvalue (schur=0, nb=0 or not included in mdrkw.m call)
%    unbalanced eigenvalue (schur=0, nb=1)
%    schur (schur=1)
% --------------------------------------------------------------------------

if (nargin<3)
   schur=1;   % Schur decomposition is  default
   nb=0;
else
   schur=comp(1);
   nb=comp(2);
end

if (track==1)
   dr=[system,'.out'];
   disp('Tracking option on')
   disp(['diary is ',dr])
   eval(['diary ',dr]);   
end

% **************************************************************************
% Step 0: Settings for computational parameters and methods
% **************************************************************************

% --------------------------------------------------------------------------
% Choose critical value for roots: mu is unstable if bcrit*mu >= 1
% --------------------------------------------------------------------------

bcrit=1/1.00001;


% **************************************************************************
% Step 1: Specifying dynamic systm                                           
% **************************************************************************

% The file lin.m, which contains the linearization and the s.s., should have
% been run already (replacing the "eval(system)" command in original
% program

eval(system)

% **************************************************************************
% Optional information on setup specified model                        
% **************************************************************************

if (track==1)
   disp('A matrix')
   disp(A)
   disp('B matrix')
   disp(B)
   disp('transformed C(F) matrix')
   %disp(CF)
   nlead=cols(CF)/nx;
   disp(['Number of leads: ',num2str(nlead)])
   disp('Predetermined endogenous variables: ')
   %disp(ylist(lpd,:))
   disp('Full vector of endogenous variables')
   %disp(ylist) 
   disp('Exogenous variables')
   %disp(xlist)
   disp('  ')
   dr=['s',system,'.mat A B CF nx lpd ylist xlist'];
  % eval(['save ',dr]); 
   disp(['Original system information saved in ',dr])
   disp(' ')
   disp('running system reduction program (redkw.m)')  
end

% **************************************************************************
% Step 2: Call to reduction program                                                
% **************************************************************************

[nd,nf,AN,BN,reord,CFN]=Redkw(A,B,CF,nx,lpd,track);

% **************************************************************************
% Optional information on output of reduction program                      *  
% **************************************************************************

if (track==1)
   disp('transformed A matrix')
   disp(AN)
   disp('transformed B matrix')
   disp(BN)
   disp('Reordering of variables')
   disp(reord)
   disp('transformed C(F) matrix')
   disp(CFN)
   disp('flow variables')
   %disp(ylist(reord(1:nf),:))
   disp('dynamic variables')
   %disp(ylist(reord(nf+1:nf+nd),:))
   disp(' ')
   disp('reordering vector')
    disp('  ')
   dr=['r',system,'.mat nd nf AN BN reord CFN']
   eval(['save ',dr]); 
   disp(['Reduced system information saved in ',dr])
   disp(' ')
   disp('running reduced system transformation program (dynkw.m)') 
end
% **************************************************************************
% Step 3: Translation matrices for fundamental dynamic equation
% **************************************************************************


[LE,MU]=Dynkw(BN,nf,ny,lpd,bcrit,track,schur,nb);

if (track==1)
   disp('reduced system translation completed')
   disp(' ')
   disp('unstable eigenvalues are: ')
   disp(diag(MU))
   disp('LE matrix that isolates unstable canonical variables is:')
   disp(LE)
   dr=['d',system,'.mat LE MU']
   eval(['save ',dr]); 
   disp(['Transformed system information saved in ',dr])
   disp(' ')
   disp('running Markov decision rule program (mdrkw.m)') 
end




% **************************************************************************
% Step 4: Solution for Markov decision rule                                        
% **************************************************************************

 

% --------------------------------------------------------------------------
%   1. Specification of forcing process                 
% --------------------------------------------------------------------------

%eval(driver)
G = eye(size(RHO));
if (0)
   disp('Driving process components (RHO,Q,G) in system specified as:')
   disp('        x(t)    = Q delta(t)')
   disp('        delta(t) = RHO * delta(t-1) + G e(t)')
   disp(RHO)
   disp(Q)
   disp(G)
end

% --------------------------------------------------------------------------
%   2. Call solution program and get output                                
% --------------------------------------------------------------------------


[M,PI]=mdrkw_learning(BN,CFN,Q,RHO,LE,MU,reord,nd,lpd,nx,bcrit,schur,bigA);

if(0)
display(['Size of M: ' num2str(rows(M)) ' by ' num2str(cols(M))])
display(['Size of PI: ' num2str(rows(PI)) ' by ' num2str(cols(PI))])
end


% --------------------------------------------------------------------------
% Augment G matrix for larger system (including endogenous variables)
% --------------------------------------------------------------------------

nlpd=max(rows(lpd),cols(lpd)); % number of endogenous states
ninv=cols(G); % number of innovations

G=[zeros(nlpd,ninv)
   G];

% --------------------------------------------------------------------------
%   Optional information on output of Markov decision rules                 
% --------------------------------------------------------------------------

if (track==1)
  disp('Markov decision rules computed')
  disp(' ')
  disp(ylist(1:ny,:))
  disp('  ')
  disp('PI matrix: y(t) elements')
  disp('  ')
  disp(PI(1:ny,:))
  disp(xlist(1:nx,:))
  disp('  ')
  disp('PI matrix: x(t) elements')
  disp('  ')
  disp(PI(ny+1:ny+nx,:))
  disp('elements of k(t)')
  disp(ylist(lpd,:))
  disp('  ')
  disp('M matrix:  elements ordered k(t), then delta(t)')
  disp('  ')
  disp(M)
  dr=['m',system,'.mat PI M G']
  eval(['save ',dr]); 
  disp(['Decision rule information saved in ',dr])
  disp(' ')
  disp('model solution completed by resolkw.m')
end

diary off






